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We develop a mean-field theory for Bose-Einstein condensation of spin-1 atoms with inter- 
nal degrees of freedom. It is applicable to nonuniform systems at finite temperatures with a 
plausible feature of satisfying the Hugenholtz-Pines theorem and various conservation laws si- 
multaneously. Using it, we clarify thermodynamic properties and the excitation spectra of a 
uniform gas. The condensate is confirmed to remain in the same internal state from T — Q up 
to T c for both antiferromagnetic and ferromagnetic interactions. The excitation spectra of the 
antiferromagnetic (ferromagnetic) interaction are found to have only a single gapless mode, 
contrary to the prediction of the Bogoliubov theory where three (two) of them are gapless. 
We present a detailed discussion on those single-particle excitations in connection with the 
collective excitations. 

KEYWORDS: Bose-Einstein condensation, mean-field theory, spinor BEC, Luttinger-Ward functional, 
Bogoliubov-de Gennes equation 



1. Introduction 

In 1998, Stamper-Kurn et al. 1 have successfully con- 
fined 23 Na atoms in an optical dipole trap. They thereby 
realized Bosc-Einstcin condensation (BEC) with inter- 
nal degrees of freedom corresponding to the three hyper- 
fine spin states \F = 1, mp = ±1, 0). Up to now, multi- 
component BEC has been observed also in other systems 
such as 23 Na (F = 2) 2 and 87 Rb (i^l 3 - 4 and J F = 2 5 ' 6 ). 
The internal degrees of freedom are expected to bring 
new physics into BEC absent in magnetically trapped 
single-component systems. 

Several theories have been constructed for the weakly 
interacting Bose gases of the single-component system. 
The basic Bogoliubov theory 7 ' 8 is a perturbation theory 
without self-consistency which is applicable only near 
T = 0. The Gross-Pitaevskii (GP) equation 9 " 11 for the 
condensate may be regarded as the inhomogeneous Bo- 
goliubov theory without the quasiparticlc contribution. 
Among extensions of the Bogoliubov theory to finite 
temperatures is the mean-field Hartree-Fock-Bogoliubov 
(HFB) theory derived with the Wick decomposition pro- 
cedure. However, it predicts unphysical energy gap in the 
excitation spectrum in contradiction to the Hugenholtz- 
Pines theorem. 12 To remove it, one generally introduces 
another approximation called the Shohno ( "Popov" ) ap- 



proximation 



13, 14 



of neglecting the anomalous pair cor- 



relation completely. However, the approximation has a 
fundamental flaw of yielding dynamical equations of mo- 
tion which does not satisfy various conservation laws. 
Thus, there are still no systematic self-consistent approx- 
imation schemes in condensed Bose systems which carry 
the two properties of the exact theory simultaneously: 15 
"gapless" and "conserving." 

Theoretical studies of multi-component BEC was 
started by Ohmi and Machida 16 and Ho 17 independently 
in 1998. They were based on the GP equation for the 
spin-1 Bose gas to determine the structure of the con- 
densate wave function at T = 0. They have also clarified 



the excitation spectra by considering fluctuations around 
the solution of the GP equation. These pioneering works 
have been followed by detailed studies on the equilibrium 
and dynamical properties of spin-1 BEC. 18-24 However, 
most of them consider only the region near T = based 
on the Bogoliubov or GP equations. In addition, other 
few studies at finite temperatures adopt the Shohno ap- 
proximation which may not provide reliable predictions 
with the reasons mentioned above, especially for multi- 
component systems. Thus, spin-1 BEC at finite temper- 
atures is theoretically not well understood yet. 

Recently, one of the authors has constructed a mean- 
field theory for a single-component BEC at finite tem- 
peratures 14, 25 which is both "gapless" and "conserving" 
with finite anomalous pair correlation. We here extend 
the theory to the three-component BEC. We then apply 
it to a uniform gas with a contact interaction to reveal 
its basic features as a function of temperature. As men- 
tioned above, no detailed studies have been made at finite 
temperatures even for the uniform system. Especially, we 
wish to clarify the condensate wave function and the ex- 
citation spectra of the multi-component system at finite 
temperatures which may have non-trivial structures. 

This paper is organized as follows. Section 2 de- 
rives a closed set of equations to determine the ther- 
modynamic equilibrium of general nonuniform three- 
component BEC. Section 3 applies the formulation to 
a uniform gas under constant density. Section 4 presents 
numerical results for spin-dependent antifcrro- and fer- 
romagnetic interaction. Section 5 concludes the paper. 
In Appendix, we present the multi-component version of 
the Hugenholtz-Pines theorem. 12 We put h = k-Q = 1 
throughout. 
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2. Formulation 

2. 1 Hamiltonian 

We consider a system described by the Hamilto- 
nian: 16, 17 

H = E/ drV£W(ffo - M)Vv(r) + \ j dr ^ 

X Wa^rr'^-^^U^^T^^r'^VV^l) ■ 

(1) 
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Here Ho = — + V cxt (r) with M the mass and V cxt 
the external potential, fj. is the chemical potential, and 
t/jaif) is the field operator with the subscript a = 0, ±1 
denoting the spin component. The quantity U aa ^ TT > is a 
spin-dependent interaction given by 

Uaa',TT'(ri - r 2 ) 

= <5(ri - r 2 ) (c <W<5tt' + c 2 F aa , ■ F TT ,) , (2) 

where Co and c 2 are parameters defined in terms of spin 
dependent s-wave scattering length a/ (/ = 0, 1, 

4tt (a + 2a 2 ) 



co 



C2 



M 3 
4tt (a 2 - a ) 



M 3 

and ;F = (-FW') is the spin matrix for F = 1: 
"0 1 0" 



F z = 



1 

V2 

/ 

71 



1 o 1 
1 



-1 



1 



1 

0-1 



■ ) as 17 
(3a) 

(3b) 

(4a) 
(4b) 
(4c) 



The second term in the round bracket of eq. (2) is sim- 
ilar in form with the Hciscnberg spin Hamiltonian. We 
hence call the c 2 > (c 2 < 0) case as antiferromagnctic 
(ferromagnetic) . 

It is convenient to introduce the Nambu vector: 



(5) 



Using eq. (5) and the normal-ordering operator A^, 26 the 
interaction in eq. (1) is transformed into 



3 

H'=^E / drdr ' u Ar - r') 



x N[^\r)A v ^[r)^{r')A v ^{r')\ . (6) 
Here IA V is defined by 

co5(r — r') : v = 



U„{r-r') = 



c 2 5(r 



v = 1,2,3 



(7) 



and A v denotes the vertex: 



A" 









(AT 

with A" defined in terms of the unit matrix 1 and eq. (4) 

as 

1 : v = 
F" :v= 1,2,3 



A v - 



(8) 



(9) 



The equivalence between the interaction term in eq. (1) 
and eq. (6) can be checked easily by writing the latter 
without using M . The expression (6) has an advantage 
that the perturbation expansion with respect to TC can 
be carried out directly by using the Nambu Green's func- 
tion; see Appendix A of ref. 14 for details. 

2.2 Mean-field equations 

To describe the condensed phase, let us express the 
field operator ip(r)= [tpi(r) ipo(r) ?/>_i(r)] T as a sum of 
the condensate wave function *&(r) and the quasiparticle 
field <j>(r) as 

i/>(r) = *(r) + <£(r). (10) 

We next introduce the Matsubara Green's function in 
the Nambu space as 

G(r,r',T) = -f 3 (T T 0(r,T)0t(r')) (11a) 



G{ry,r) 
-F*(r,r',r) 



£(r,r',r) 
-G*{r,r',r) 



(lib) 



with T3 denoting the third Pauli matrix in the Nambu 
space. Equation (lib) can be Fourier-transformed as 



6(r,r',r) = jjj^r.r', 



iuj n )e~ 



(12) 



where (3 = T 1 and uj n is the Matsubara frequency uj n = 
2mr/f3 with n = 0, ±1, ±2, ■ • • . 

Using eq. (12), we now write down our mean-field 
Luttinger-Ward functional for BEC as 14 ' 27 



n = J dr&(r)(H - 



+ ^ ^ Tr \n(K + £ - icjJ) + G± 



(13) 

Here K is defined by K = f 3 (H — fi), i(iw„) denotes 



i(w n ) 



le^" 0+ 







le 



-iu;„Cf 



(14) 



with + an infinitesimal positive constant, and Tr in- 
cludes integration over space variables. The quantity S 
is the proper self-energy obtained from the functional 
$ = G] by 

E(r, r', u n ) = -2/3-,— -. (15) 

5G(r',r,L0 n ) 

With eq. (15), functional fl is stationary with respect to 
a variation of G which satisfies Dyson's equation: 



CT 1 = iuj n l - K 



(16) 
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The condensate wave function \& in equilibrium is also 
determined by Sfl/S^* = 0. Noting SQ/SG = 0, we ob- 
tain 

5$ 



(Ho - 



5**(r) 



0. 



(17) 



Equations (15) and (17) constitute the ^-derivative ap- 
proximation where various conservation laws are obeyed 
automatically. This is one of the main advantages of us- 
ing the Luttinger-Ward functional. 

Now, we choose $ so that the Hugcnholtz-Pines the- 
orem is satisfied simultaneously; see Appendix for the 
multi-component version of the Hugcnholtz-Pines theo- 
rem. Explicitly, our $ is given by 

$ = - E / drdr'U v (r - r') 

2 v J 

x {^Tr[A v 4>(r)4?\r)]^Tr[A v 4?(r')&(r')} 



x -Tr[A y f 3 G(r',r',z Wn )lK) 



-Tr[A v f 3 if(r)^(r')A v G(r',r,cj n )i(u; n )] 
l]r(I T r[i^ 3 G(r,r,^)] 



x -Tr{A»T 3 G(r',r',u m )} 



+ 



^Tr[A»G(r, r', w n )i(w n )A v G(r' , r, w m )i(w m )]) }. 



(18) 



This functional <f> is different from $hfb of the HFB 
theory in the Fock terms, i.e., the third and the fifth 
terms on the right-hand side of cq. (18). Indeed, ^hfb 
is recovered from $ by replacing f 3 \&\I't — > and 
G — > f 3 G. Put it another way, $ is obtained from $hfb 
by the following subtraction: 

$ =<I>hfb ~ E / drdr'W 1/ (r - r') 

x ( - ^ E Tr[i^(r)*t (r ')i^ 3 F(r', r, w n )i(w n )] 
n 

+ iE Ti '[ i ^( r ' r >») i (^) 

x i^f 3 F(r / ,r,cj m )i(a; m )] 



with 

F(r,r',w n ) = 



£(r,r',o; n y 
-F*(r,r>„) 



(19) 
(20) 



This subtraction has the effect to cancel the overcounting 
of the anomalous pair correlation in the HFB theory, 
thereby leading to a mean-field theory which satisfies the 
Hugenholtz-Pines theorem. 



Now, the proper self-energy is obtained with eqs. (15) 
and (18) as 

t(r r')-\ Ury) A(r ' r,) ' 
1 ' ~ [- A* (r , r') -£* (r , r') 

= E{^( r - r ')^ f3 / dr "UAr ~ r") 

x Tr[A»t 3 p(r",r")] + U v (r - r')A u p(r,r')A l 



(21) 

where p is the one-particle density matrix in the Nambu 
space defined by 



p(r,r') 

= f 3 #(r)* t (r' 



p(r,r') p(r,r') 
-p*(r,r') -p*(r,r') 



(qp, 



(r,r') 



with 



) (qp) (r,r') 



^G(r,r', Wn )lK). 







(22) 
(23) 



The submatrices £ and A satisfy S*(r,r') = £ (r',r) 
and A(r, r') = A T (r', r). Similarly, we have p*(r, r') = 
p T (r',r) and p T (r,r') = p(r',r). 

Next, eqs. (17) and (18) lead to the the equation for 
the condensate wave function: 

(H - M )*(r) 

+ J dr%(r,r')V(r) - A(r,r')**(r')] = 0. (24) 

Equation (24) is the GP equation in our mean-field the- 
ory. In the homogeneous case, we can write iff— > y / n u T7, 
where n is the condensate density and rj (|tj| = 1) spec- 
ifies the direction of the order parameter. Equation (24) 
then reduces to pC — £fc=oC — A fc=0 £*, which is the 
Hugenholtz-Pines relation of the multi-component case; 
see eq. (A-9) in Appendix. It tells us the existence of 
a single gapless excitation spectrum in the system, as 
shown explicitly below eq. (39). Note that p and r\ should 
be determined here as the lowest eigenvalue and its eigen- 
state, respectively. Finally, the thermodynamic relation 
N = —dil/dp yields the expression for the total particle 
number as iV = Tr J p(r,r)dr. 

In order to diagonalizc Green's function which satisfies 
eq. (16), we consider the Bogoliubov-de Gcnnes (BdG) 
equation: 



dr'H(r,r') 



u 3 (r') 
-v*(r>) 



-v*(r) 



(25) 



where H(r,r') = K(r)5(r — r') + £(r,r'). We assume 
that the eigenstate for Ej > can be normalized as 

(26) 



= s 



Jj' ■ 



It then follows from the symmetry of H that the eigen- 
vector with the eigenvalue — Ej is given by [— vj , it^] T . 14 

Let us introduce matrices U and E by 

"J , (27) 



" u 


-V 


, E = 


~E " 


-V* 


u* 


-E 
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with U=0i 



■),V= {v jl v h ■••) and (E) 3 



EjSjji. Using these matrices, the BdG equation can be 
written compactly as / dr'H(r,r')U(r') = U(r)E. It 
hence follows that we can expand Green's function as 



G(r,r',uj n ) = UirHiuJ-fy-Wir'). 



(28) 



We now carry out the summation over the Matsubara 
frequency in eq. (23). We then obtain the quasiparticle 
contribution to the density matrix in eq. (22) as 

£ («Jp)( P ,r , ) = (T T 0(r,O-)^(r')> 

= {mr)f(muHr') + V(r)[l+f(m}VHr')}, (29a) 

p^(r,r / ) = (0(r)0(r / )) 
= l{v(r)[l+ 2f(E)] LL T (r') 

+ U(r)[L+2f(E)}V T (r')}, (29b) 

where / is the Bose distribution function. 

Equations (24) and (25) with eqs. (21), (22) and (29) 
form a closed set of self-consistent equations to determine 
the equilibrium. 

Using eqs. (24) and (25), we can transform the ther- 
modynamic potential of eq. (13) into 

« cq = /r 1 J2 Mi - e"^) - E^7 M r )l 2 dr 

3 3 

- i J drdr'Tr [E(r, r')p(r', r) - A(r, r')p(r', r)] . 

(30) 

Also, entropy is obtained as 

s = - E ln /(^) - 1 1 + mi + f( E i)]} ■ 

(31) 

Finally, the specific heat is calculated by the thermody- 
namic relation C = -TdS/dT. 

3. Uniform Gas under Constant Density 

3.1 Equations to determine equilibrium 

We now apply the previous formulation to a uniform 
Bose gas with fixed density n = N/V, where V denotes 
volume of the system. In this case, the condensate wave 
function \& can be written as 



nor), 



(32) 



with no = No/V the condensate density and r\ a constant 
vector with \r\\ = 1. Also, eq. (22) may be expanded in 
plane waves as 



p(r,r') 
with 



1 E 

v 4^ 



p(k) 
-p(-k) 



p(k) 
-p(-k) 



p(fc)=WoW + p ,,p) W, 

p(k) = 5 ko N r]r] T 



p (qp) (fe). 



(33) 

(34a) 
(34b) 



The transition temperature of the nonintcracting three- 
component BEC is given by 



In 



2tt 
M 



3C(|; 



(35) 



where £(§) ~ 2.612 is the Ricmann ^-function. This To 
is 3 -2 ^ 3 times the transition temperature of the one- 
component BEC. 

The self-energy (21) becomes fc-independent in mo- 
mentum space for the contact interaction of eq. (7) . In- 
deed, E and A in k space are obtained as 

2/3 

E S v [A v Tvp(k) A" + A v p{k)A v ] , 



2 

N 



'I 5 



(36a) 



A 



2 

N 



3( 



2/3 



E^"£(*)U") T . (36b) 



where b v denotes dimensionless interaction strength de- 
fined in terms of eq. (3) by 



5u 



M , 

— c n 1/3 : v = 

47T 



M 

47T 



(37) 



c 2 n 



1/3 



i/ = 1,2,3 



The GP equation (24) is transformed by using eq. (32) 
into 

m = St/- At?*. (38) 

This is the Hugenholtz-Pines relation of the three- 
component system; see Appendix for details. The BdG 
equation (25) now reads 



Hi 



—v" 



ka 



= E, 



kcr 



Uk 

—v' 



(39) 



where H k — TS ~ AH + E, the eigenstate is normal- 



ized as \uko 



\Vko 



1, and the subscript a = 0, ±1 



distinguishes the three eigenvalues of each k. This eigen- 
value problem includes eq. (38) as the special case for 
-E-fc=OCT = with itfc=OCT = i'fc=OCT = C- It hence follows that 
the quasiparticle excitation is gapless as k — > at least 
for one branch among <r = 0, ±1. The quasiparticle con- 
tribution (29) to the density matrix is now given by 

£ foP)(fc) = {u(k)f(E k )U\k) + V(k)[l+f(E k )]v!(k)} , 

(40a) 



p^(k) = -{V(k)\l + mmiL T (k) 

+U(k)[l+2f(E k )}V T (k)}, (40b) 



with U(k)=[uki u k0 Uk-i] and V(k) = [v k i v k0 Vk-i]- 
Finally, wc sum eq. (34a) over k and take trace of it 
subsequently. We thereby obtain an expression for the 
condensate density as 



"0 



1 



(qp) 



(k). 



(41) 
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3.2 Numerical procedures 

Numerical calculations of the above self-consistent 
equations have been carried out as follows. 

(i) The chemical potential [i and the vector rj are de- 
termined by eq. (38) as the smallest eigenvalue and its 
eigenstate. To carry this out, we have to know the ex- 
pressions of no, p( qp \k) and /r qp '(fe) in advance, as seen 
from eqs. (34) and (36). We use the values from the pre- 
vious calculation for them; we initially set no = n and 
p(qp) =p(iP) =o to start the whole calculation at T = 0. 
We determine r\ self-consistently so that it actually cor- 
responds to the eigenstate of the smallest eigenvalue. 

(ii) Using /i and £ thus calculated, we solve the BdG 
equation (39) to obtain the quasi-particle spectrum E ka 
and the matrix 

(hi) With E ka and U k , we compute £>( qp) (fc), /5 (qp) (fc) 
and no by eqs. (40a), (40b) and (41), respectively. 

We iterate procedure (i)-(iii) until a convergence is 
reached. We start the whole calculation from T = and 
increase the temperature gradually. We identify T c nu- 
merically as the point where we can no longer find a 
solution of positive no. 

When we calculate A for T < T c by using eqs. (36b) 
and (40b), we observe that the summation of p^ qp ^(fe) 
over k diverges. In order to remove it numerically, we 
have introduced the energy cutoff e c so that it satisfies 

l<e c <0.3^ 2 . 14 

For T > T C) we put n = 0, A = and F = 0. The 
calculation procedures (i)-(iii) are nearly the same except 
that the chemical potential is now determined by 

Ar = E e xp[^ fcCT ( u )]-l- (42) 

4. Numerical Results 

4-1 Thermodynamic properties 

We have fixed the interaction parameters of eq. (37) as 
5q = 0.0075 and C2/co = 0, ±0.1. We have confirmed that 
the stable state is given by r?= [010] T (rj= [100] T ) at all 
temperatures for the antiferromagnetic (ferromagnetic) 
interaction of C2/co = 0.1 {c^/cq = — 0.1). We have found 
for all these cases that the transition temperature T c is 
increased over To by (T c — Tq")/Tq ~28o, which is of the 
same order as (T c — To) /To — 2.335 of the single-component 
case. 14,25 We have also observed that T c increases mono- 
tonically as a function of C2/C0. 

Figure 1 shows temperature dependence of the con- 
densate density in comparison with the ideal Bose gas 
result no/n = 1 — (T/T c ) 3 / 2 . We observe that no becomes 
larger as the interaction is increased. We have also con- 
firmed that n depends little on C2 as long as |ca/co| is 
much smaller than 1. 

Figure 2 displays temperature dependence of the spe- 
cific heat per particle in comparison with the ideal gas 
result. As seen clearly, a divergent behavior shows up 
just below T c for the finite interaction. It is quite similar 
to the behavior of the single-component system. 

Figure 3 plots temperature dependence of the chemical 
potential scaled by To. With finite interaction, we find a 
small peak with a discontinuity at T = T C . This disconti- 
nuity is a signal of a first-order transition and also seen 




0.2 0.4 0.6 0.8 1 
T/T c 



Fig. 1. Temperature dependence of density of particle. 




T/T c 



Fig. 2. Temperature dependence of specific heat. 




0.02 1 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 1.2 

T/T c 

Fig. 3. chemical potential 

in the single-component case. 

4-2 Energy dispersion of antiferromagnetic interaction 
We next study the energy dispersion of the antifer- 
romagnetic interaction ic-ijco = 0.1) at T = and com- 
pare our result with the prediction of the Bogoliubov ap- 
proximation (BA). In BA, the quantities /1, r\ and E ka 
can all be obtained analytically as /i = cono, t/= [010] T , 
Eko = V^feI^fe+2cono) and E k ±i = V £ fc( £ fc + 2c 2 ri o) with 
e/j = fc 2 /2M; thus, Ek±i is degenerate in BA. 

Figure 4 shows energy dispersion at T = 0. The BA 
predicts that all the three branches have /c-lincar sound- 
wave dispersions in the long wavelength limit and two 
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Fig. 4. The energy dispersion of antiferromagnetic interaction at 
T = 
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Fig. 5. Temperature dependence of excitation energy gap in anti- 
ferro interaction. 

of them are degenerate. Our theory also yields a sin- 
gle gapless excitation in agreement with the Hugcnholtz- 
Pines theorem, and its sound velocity caf coincides al- 
most completely with the expression cquq /M of B A for 
the non-degenerate dispersion, as seen in Fig. 4. As for 
the other two branches, our mean-field theory predicts 
degenerate gapful excitations. This gap is brought about 
by the anomalous pair potential of the quasiparticlc field 
in the BdG equation, which is absent in BA. 

The qualitative features of the energy dispersions men- 
tioned above remain invariant up to T c . The sound ve- 
locity decreases monotonically towards T c . In contrast, 
the energy gap between -Efc = o,o and -Efc=o,±i displays a 
more complicated behavior. Figure 5 shows temperature 
dependence of the energy gap. The gap has a local maxi- 
mum as a function of temperature, which may be traced 
to the temperature dependence of the anomalous pair 
correlation ^2 k p(fc)( qp ). Indeed, ^2 k /5(fc)( qp ' has a local 
maximum even in the single-component case, which is 
due to the competition between decreasing A and in- 
creasing Bosc distribution function as T is raised. 

We finally comment on the finite energy gap of 
i?fc = o,±i predicted by our mean-field theory in con- 
nection with the collective excitations. Gavoret and 
Nozicrcs 28 performed a detailed study on the structure of 
the perturbation expansion for the one-component BEC 
at T = 0. They have thereby established that the single- 
particle (i.e., Bogoliubov) and collective excitations share 



common dispersions. Indeed, the statement holds even at 
any finite temperatures of T <T C , which stems from the 
fact that the single-particle propagator appears in the in- 
termediate state of the collective propagator in the con- 
densate; see eq. (2.33) of Szepfalusy and Kondor 29 and 
also their discussion in the paragraph below eq. (2.34). 
However, it should be noted that this theorem only states 
that the single-particle excitation spectra are also in- 
cluded in the collective excitations, but not vice versa. 
Thus, it is neither surprising nor contradictory that the 
collective excitations have extra branches not present in 
the single-particle excitations. 

The theorem may be extended to multi-component 
systems to tell us that all the single-particle excitation 
spectra are also present in the collective excitations. It 
is also worth noting that the Hugenholtz-Pines theorem 
of the multi-component case only states the existence 
of a single gapless one-particle excitation spectrum, as 
shown in Appendix. Thus, the qualitative features of the 
three branches E k0 and E k ±i obtained here satisfy all 
the exact statements, and it is not contradictory that 
the collective spin-wave modes are not included in the 
single-particle excitations. 

Recently, Szepfalusy and Szirmai 30 have extended the 
consideration of Szepfalusy and Kondor 29 to the spin- 
1 system. They conclude that there appear additional 
collective excitations due to the spin-degrees of freedom 
whose spectra in the long-wavelength limit are also gap- 
less and identical with those of the corresponding single- 
particle excitations. This result is clearly in contradiction 
with the finite energy gap of -Efc=o,±i in our mean-field 
theory. We suspect that the conclusion of Szepfalusy and 
Szirmai 30 on the existence of the spin-wave mode in the 
single-particle spectra is due to an inappropriate assump- 
tion on the poles and zeros of the single-particle Green's 
function expressed with the same denominator as the 
collective propagators. Note also that they only confirm 
their general theory within the Bogoliubov approxima- 
tion and the random-phase approximation of retaining 
only the Hartree term where the anomalous pair poten- 
tial due to quasiparticles was neglected completely. 

4-3 Energy dispersion of ferromagnetic interaction 

We move on to study the energy dispersion of the fer- 
romagnetic interaction (c^Icq = — 0.1). In this case, BA 
near T = predicts = (c + c 2 )n , r) = [100] T , E k+ i = 
v /e fc [e fc + 2(co + c 2 )no]- E k0 = e k and E k -i = e k -2c 2 n . 

Figure 6 shows the energy dispersions of E k+ \ and 
E k o at T = 0. According to BA, the two branches are 
both gapless and E k+1 (E k0 ) has a linear (quadratic) 
dispersion in the long wavelength limit. However, our 
mean-field theory predicts a gapless excitation only for 
.Efc+i, and the other branch E k g acquires an energy 
gap. The existence of a gapless excitation E k+ i agrees 
with the Hugenholtz-Pines theorem, and its sound veloc- 
ity cfe coincides almost completely with the expression 
(cq + C2)no/M of BA. On the other hand, a finite gap in 
E k o is brought about by the anomalous pair potential of 
the quasiparticlc field. 

We have confirmed that the vector r] remains invariant 
at finite temperatures. Also, the qualitative features of 
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Fig. 6. Energy dispersion of ferromagnetic interaction at T = 
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Fig. 7. Temperature dependence of excitation energy gap in ferro 
interaction. 



the dispersions do not change up to T — T c . The sound 
velocity decreases monotonically towards T c . Figure 7 
shows the energy gaps of Ej~q and -Efc-i at k = as a 
function of temperature. The gap of Eko has a maximum 
at finite temperature, which has the same origin as in the 
antiferromagnetic case. 

5. Summary 

We have extended a conserving gapless mean- field the- 
ory of single-component Bose-Einstein condensates to a 
three-component system. The equations to determine the 
equilibrium are given by eqs. (24) and (25) with eqs. (21), 
(22) and (29), which are applicable to nonuniform sys- 
tems such as the trapped BEC with and without vortices. 
They arc compatible with the Hugcnholtz-Pincs theorem 
and may be extended easily to general multi-component 
systems. Another advantage of our theory is that various 
conservation laws are automatically obeyed when applied 
to dynamical phenomena. Thus, the static and dynamical 
properties of BEC may be treated systematically within 
a single theoretical framework. 

We have applied the mean-field theory to a uniform gas 
with constant density. A qualitative difference is found 
in the quasiparticle dispersions between our theory and 
the Bogoliubov theory near T = 0. The Bogoliubov ap- 
proximation predicts that there are three (two) gapless 
excitations for the antiferromagnetic (ferromagnetic) in- 
teraction. On the other hand, our theory yields only a 



single gapless branch with a sound-wave dispersion at 
low energy, and the others are pushed up to higher en- 
ergies due to the self-consistent potential of the quasi- 
particle field. We have presented at the end of §4.2 a 
detailed discussion on those single-particle excitations in 
connection with the collective excitations. The speed of 
sound obtained here well agrees with the corresponding 
branch in the Bogoliubov approximation. We have also 
clarified finite-temperature effects on the quasiparticle 
dispersions as well as temperature dependences of basic 
thermodynamic quantities. 
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Appendix: Hugenholtz-Pines theorem for spin-1 
BEC 

We here generalize the Hugenholtz-Pines theorem 12 to 
the spin-1 BEC by closely following the procedure of 
Popov. 31 It will be shown that there is at least a sin- 
gle gapless excitation spectrum in the spin-1 BEC (and 
also in any multi-component BEC). 

First of all, we present Feynman rules for the bare 
perturbation expansion of the spin-1 condensate as the 
basis for our proof. They are obtained by slightly mod- 
ifying the rules of the single-component system. To see 
this, let us rewrite the interaction of eq. (6) as 

3 

W = ^ f drdr'U u {r - r') 

x Af[^(r)A v ip(r)^(r')A v ip(r')} , (A-l) 

where A" is given by eq. (9), and the field operators are 
defined by 



(A-2) 



The expression (A-l) enables us to perform the per- 
turbation expansion of the spin-1 system in terms of 
the Feynman diagrams of the single-component sys- 
tem, i.e., no new diagrams are necessary. To be more 
specific, we first consider the normal system with- 
out the condensation. We then notice that: (i) the 
pairs ij)^ (r)A^ xp(r) and ip' (r')A v ip(r') can be moved 
around anywhere within the Af and/or T T operators 
in the perturbation expansion; and (ii) a contraction 
of (i = T*jTi) with V^(j) automatically yields the 

3x3 matrix (T T ip(i)^ (j))o = — G (i,j), where the 
subscript denotes the average with respect to the 
non-interacting Hamiltonian. Also, the final contraction 
within a closed particle loop can be transformed as 
(T T f(f)M(i,j)1>V))o = TrM(i,3)(T T xP(j)^(i)) with 
M_(i,j) denoting some matrix product of contractions. 
We now realize that the expansion can be performed in 
terms of the Feynman diagrams of the single-component 
system with the following modifications: (i) Associate 
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G (i,j) for the particle line from j to i and the ma- 
trix A v for each vertex, and (ii) take Tr for every closed 
particle line and perform summations over v. 

We next consider the condensed system and adopt the 
shift transformation (10) with the self-consistency con- 
dition: 



a 2 n 



<<Kr)) = 0. 



(A3) 



Let us substitute eq. (10) into eq. (1) and carry out the 
perturbation expansion by choosing 



(A4) 



as the free field. Once again, we only need the Feynman 
diagrams of the single-component system, and the Feyn- 
man rules for the scalar quantities Go, ^ and m are 
modified as follows: (i) G (i,j)— >G (i, j), *(r) —> *(r), 
and m(r) — ► **(r); (ii) A" for the vertex; (iii) take Tr 
for every closed particle line and perform summations 
over v. 

Now, we can directly extend the consideration of 
Popov 31 on the general structure of the perturbation ex- 
pansion to the spin-1 system. We thereby conclude that 
eq. (A-3) is equivalent to the extremal property: 

(5**(r) v ' 

It is also straightforward to prove the following relation: 



s 2 n 



8* a .(r')6**(r) 



>(r,r';w n = 0) - fi5 aa <5(r - r'), 

(A- 6a) 



s 2 n 



= A CTCT /(r,r';cj„ = 0). 



(A-6b) 



J**,(r')<^*(r) 

We finally consider the uniform condensate, write the 
condensate wave function as eq. (32), and expand S and 
A as eq. (33). It is then evident that the thermodynamic 
potential depends on & and via the product tiq = 
\ft\I>. Hence we can write eqs. (A-5) and (A-6) for the 
uniform case as 

asi on 

dm~ "cW 



(A7) 



and 



d 2 n 



on 



d 2 n 



d^^dm 



= <W ^— + K'^<t^, ( A ' 8a ) 



dm, am 



d 2 n 

dnl '' 



(A-8b) 



respectively. Using eqs. (A-5)-(A-8), we arrive at the gen- 
eralized Hugenholtz-Pines theorem: 

[E <T<T> (fc = 0,w n = 0)* <T / - A aa ,{k = 0,cj n = 0)m a ,] 



E 



(A-9) 



It tells us that there is at least a single gapless one- 
particle excitation in the spin-1 BEC; see the discussion 
below eq. (39) on this point. It is also straightforward 

to extend the consideration to any multi-component sys- 
tem, which leads to the same conclusion. 
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